Coral health statuses were monitored for disease in the CCMI nursery Dataset is located in Dryad: Brown, Anya et al. (2022), CCMI nursery coral disease 2019, Dryad, Dataset, https://doi.org/10.25338/B8F643
library(tidyr)
library(dplyr)
library(reshape2)
library(plyr)
library(ggplot2)
library(Rmisc)
library(stringr)
library(lme4)
library(lmerTest)
library(cowplot)
library(patchwork)
library(viridis)
library(glmmTMB)
library(DHARMa)
#dataset in Dryad
Data set up
dis_1 <-read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Frame Surveys/Data/FrameSurvey_Dec2019_wideformat_edit2.csv")
setwd("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Frame Surveys/Data/")
#if goes from dead to missing, change back to dead
library(data.table)
Make sure data in columns are consistently named
dis <- reshape2::melt(dis_1, id= c("Frame", "Column","Row","Genotype")) #value.name ="Coral_stat")
colnames(dis)[5] <- "Date"
colnames(dis)[6] <- "Coral_stat"
#remove X from the dates
dis$Date <-gsub("X","",dis$Date)
#create unique IDs for each coral based on frame number and location on a frame
dis$ID <- paste(dis$Frame, dis$Column, dis$Row)
#remove Xs from data
#dis$Date <- str_remove(dis$Date,"X")
#remove leading or trailing white spaces
dis$Coral_stat <- trimws(dis$Coral_stat)
#Replace P with H
p <- which(dis$Coral_stat == "P")
dis[p,"Coral_stat"] <- "H"
#Replace XS with X
s <- which(dis$Coral_stat == "XS")
dis[s,"Coral_stat"] <- "X"
#Replace HR with H (Healty recovred should be just health)
h <- which(dis$Coral_stat == "HR")
dis[h,"Coral_stat"] <- "H"
#replace "0" with O for missing
o <- which(dis$Coral_stat == "0")
dis[o,"Coral_stat"] <- "O"
#replace OD with Dead
d <- which(dis$Coral_stat == "OD")
dis[d,"Coral_stat"] <- "D"
#replace R for Recovered with Healthy
e <- which(dis$Coral_stat == "R")
dis[e,"Coral_stat"] <- "H"
#replace diseased and recovered (XR) with X (disease)
f <- which(dis$Coral_stat == "XR")
dis[f,"Coral_stat"] <- "X"
#replace recovered with disease with just disease
g <- which(dis$Coral_stat == "RX")
dis[g,"Coral_stat"] <- "X"
#remove data from May 16 - not trustworthy
dis <- dis[which(dis$Date != "2019.05.16"),]
#remove missing corals (never present) and where there is no data
dis <- dis[which(dis$Coral_stat != "O" & dis$Coral_stat != "NoData"),]
Genotype counts
counts_dis_g <- plyr::ddply(na.omit(dis), c("Date", "Genotype","Frame","Coral_stat"), summarize,
count = length(Coral_stat))
Including zeroes for other health categories
#put all data longwise
counts_dis_g <- as.data.table(counts_dis_g)
dim(counts_dis_g)
## [1] 2668 5
counts_cast <- dcast.data.table(counts_dis_g, Date+Genotype+Frame~Coral_stat, value.var = "count")
counts_cast <- counts_cast %>% mutate_all(~replace(., is.na(.), 0))
dim(counts_cast)
## [1] 1433 6
Melt data
counts_cast_melt <- melt(data = counts_cast, id.vars = c("Date","Frame","Genotype"))
colnames(counts_cast_melt)[4] <- "Coral_stat"
colnames(counts_cast_melt)[5] <- "count"
Labeling single and mixed genotype frames
library(data.table)
frame <- data.table(unique(counts_cast_melt$Frame))
colnames(frame)<- "Frame"
frame$sing_mix <- ifelse(frame$Frame > 17, "mixed","single")
s <- which(frame$Frame == "1207")
frame[s,"sing_mix"] <- "single"
w <- which(frame$Frame == "1279")
frame[w,"sing_mix"] <- "single"
Counts across all of the frames and genotypes
#count number of corals in a health category - this sums the number of corals over all health categories on a frame to get the total number of corals on a frame
counts_fr <- ddply(counts_cast_melt, c("Date", "Frame"), summarize, countf = sum(count))
#for each frame, wheter it's a single or mixed genotype
singmixf <- merge(frame, counts_fr)
#total number of corals on single or mixed frames at each date
total.per <- ddply(singmixf, c("Date","sing_mix"), summarize,
total = sum(countf))
#total number of frames
total.f <- ddply(singmixf, c("Date","sing_mix"), summarize,
total = length(Frame))
#12 mixed frames
#18 single frames
#432 single corals
#218 mixed
Calculation of proportion of health statuses on a frame including zeroes
counts_all_singmix <- merge(counts_cast_melt, singmixf, by = c("Frame", "Date"))
dim(counts_all_singmix)
## [1] 4299 7
counts_all_of_frame <- ddply(counts_all_singmix, c("Date", "Frame","sing_mix","countf","Coral_stat"), summarize,
count_all = sum(count))
dim(counts_all_of_frame)
## [1] 1872 6
#average over all frames
sum_all_counts <- Rmisc::summarySE(counts_all_of_frame, measurevar = "count_all", groupvars = c("Date","Coral_stat","sing_mix"))
counts_all_of_frame$prop <- counts_all_of_frame$count_all/counts_all_of_frame$countf
#average of proportion in each health category per frame
sum_all <- Rmisc::summarySE(counts_all_of_frame, measurevar = "prop", groupvars = c("Date","Coral_stat","sing_mix"))
Figure 2: Plot of Disease for all corals across the nursery
sum_all2 <- subset(sum_all, Date != "2019.02.01")
sum_all2$Date2 <- gsub("[.]", "/", sum_all2$Date)
sum_all2$Date3 <- as.Date(sum_all2$Date2)
#### Used in Figure 2 below ####
sing_mix_all <- ggplot(sum_all2[which(sum_all2$Coral_stat=="X"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5)) + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom") + xlab("Date")
Disease-only data
disease <- counts_all_of_frame[which(counts_all_of_frame$Coral_stat == "X"),]
dim(disease)
## [1] 624 7
dead <- counts_all_of_frame[which(counts_all_of_frame$Coral_stat == "D"),]
total.disease.counts <- ddply(disease, c("Date","sing_mix"), summarize,
sum.disease = sum(count_all), sum.all = sum(countf))
total.dead.counts <- ddply(dead, c("Date","sing_mix"), summarize, sum.disease = sum(count_all), sum.all = sum(countf))
#create factors and turn dates into date data type
disease$Frame <- as.factor(disease$Frame)
disease$Date2 <- gsub("[.]", "/", disease$Date)
disease$Date3 <- as.Date(disease$Date2)
#disease$Date2 <- as.numeric(disease$Date)
#includind a polynomial
x2 <-poly(disease$Date3,2)
Model 1 analysis (all data)
disease <- counts_all_of_frame[which(counts_all_of_frame$Coral_stat == "X"),]
disease$Frame <- as.factor(disease$Frame)
disease$Date2 <- gsub("[.]", "/", disease$Date)
disease$Date3 <- as.Date(disease$Date2)
#include density
Frame_area <- read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Frame Surveys/Data/Frame_area.csv")
#dim(disease)
#dim(disease_dens)
disease_dens <- merge(disease, Frame_area, by = "Frame")
disease_dens$density <- disease_dens$countf/disease_dens$Area
disease_dens$density_rescale <- reshape::rescaler(to = c(0,1), x = disease_dens$density)
Model 1: Binomial analysis
x2 <-poly(disease_dens$Date3,2)
library(glmmTMB)
glmmtb_alldat <- glmmTMB(cbind(count_all, countf-count_all) ~ sing_mix*x2+density + (1|Frame), family = binomial, data= disease_dens)
Check Model 1 residuals
library(DHARMa, quietly = TRUE)
glmtb_resid1 <- simulateResiduals(glmmtb_alldat)
testOutliers(glmtb_resid1, type = "bootstrap")
##
## DHARMa bootstrapped outlier test
##
## data: glmtb_resid1
## outliers at both margin(s) = 0, observations = 624, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01686699
## sample estimates:
## outlier frequency (expected: 0.00259615384615385 )
## 0
plot(glmtb_resid1)
Summary of Model 1 as a table
summary(glmmtb_alldat)
## Family: binomial ( logit )
## Formula:
## cbind(count_all, countf - count_all) ~ sing_mix * x2 + density +
## (1 | Frame)
## Data: disease_dens
##
## AIC BIC logLik deviance df.resid
## 2997.4 3032.9 -1490.7 2981.4 616
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Frame (Intercept) 2.333 1.527
## Number of obs: 624, groups: Frame, 30
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6730 0.5349 -3.128 0.00176 **
## sing_mixsingle 1.0493 0.6604 1.589 0.11209
## x21 11.3317 2.1391 5.297 1.17e-07 ***
## x22 -23.3456 2.4903 -9.375 < 2e-16 ***
## density -0.1371 0.1097 -1.250 0.21140
## sing_mixsingle:x21 12.1196 2.6992 4.490 7.12e-06 ***
## sing_mixsingle:x22 -8.5804 3.1386 -2.734 0.00626 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(glmmtb_alldat))
| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| sing_mix | 2.658957 | 1 | 0.1029682 |
| x2 | 358.706944 | 2 | 0.0000000 |
| density | 1.561784 | 1 | 0.2114042 |
| sing_mix:x2 | 23.471265 | 2 | 0.0000080 |
Density information
disease_dens %>%
dplyr::group_by(Date,sing_mix)%>%
dplyr::summarize(total = sum(count_all))
## `summarise()` has grouped output by 'Date'. You can override using the
## `.groups` argument.
Figure 1 data summary
sum_all2 <- Rmisc::summarySE(counts_all_of_frame, measurevar = "prop", groupvars = c("Date","Coral_stat"))
sum_all2$Date2 <- gsub("[.]", "/", sum_all2$Date)
sum_all2$Date3 <- as.Date(sum_all2$Date2)
Figure 1 Plot
#### Figure 1 ####
ggplot(sum_all2[which(sum_all2$Date!="2019.02.01"),], aes(x = Date3, y = prop, group = Coral_stat)) + geom_point(aes(fill = Coral_stat),color = "black", pch = 21, size = 3) + theme_bw() + geom_errorbar(aes(ymax = prop+se, ymin = prop-se)) + ylab("Proportion of health status over time") + theme(axis.text.x = element_text(angle = 90), panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + scale_fill_manual(values = c("black","orange", "gray"), name = "Coral Health", labels = c("Dead","Healthy", "Disease")) + xlab("Date")
Figure S2: Health statuses across all frames set up
#### Supplement - all Frames
sum_all3 <- Rmisc::summarySE(counts_all_of_frame, measurevar = "prop", groupvars = c("Date","Coral_stat", "Frame", "sing_mix"),na.rm = T)
## Warning in qt(conf.interval/2 + 0.5, datac$N - 1): NaNs produced
sum_all3$Date2 <- gsub("[.]", "/", sum_all3$Date)
sum_all3$Date3 <- as.Date(sum_all3$Date2)
sum_all3$Frame <- factor(sum_all3$Frame, levels = c("1235","1231","14","9","1237","7","5","1221","1219","12","11","1251","1201","1217","13","10","1209","6","4","3","8","1279","1207","16","15","1107","1","2","1295","1296"))
Figure S2: Health statuses across all frames
ggplot(sum_all3[which(sum_all3$Date!="2019.02.01"),], aes(x = Date3, y = prop, group = Coral_stat)) + theme_bw() + geom_errorbar(aes(ymax = prop+se, ymin = prop-se)) + ylab("Proportion of health status over time") + theme(text = element_text(size = 20), axis.text.x = element_text(size = 16, angle = 90), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background = element_rect(color = "black", fill ="white")) + scale_color_manual(values = c("black","orange", "gray"), name = "Coral Health",labels = c("Dead","Healthy","Disease")) + xlab("Date") + facet_wrap(~Frame+sing_mix) + geom_line(aes(color = Coral_stat), size = 1)+ geom_point(aes(color = Coral_stat), size = 1)
Including genotype information calculates the proportion of disease within a genotype on a frame
counts_cast_melt <- melt(data = counts_cast, id.vars = c("Date","Frame","Genotype"))
colnames(counts_cast_melt)[4] <- "Coral_stat"
colnames(counts_cast_melt)[5] <- "count"
counts_cast_melt <- subset(counts_cast_melt, Coral_stat == "X" | Coral_stat == "H" | Coral_stat == "D")
dim(counts_cast_melt)
## [1] 4299 5
#Blue
B <- counts_cast_melt[which(counts_cast_melt$Genotype == "B"),]
counts_fr <- ddply(B, c("Date", "Frame"), summarize,
countf = sum(count))
Bf <- merge(B, counts_fr, by = c("Date","Frame"))
Bf$prop <- Bf$count/Bf$countf
genotypefunc <- function(Gen){
Gent <- counts_cast_melt[which(counts_cast_melt$Genotype == Gen),]
counts_fr <- ddply(Gent, c("Date", "Frame"), summarize,
countf = sum(count))
Gentf <- merge(Gent, counts_fr, by = c("Date","Frame"))
Gentf$prop <- Gentf$count/Gentf$countf
Gentf <- as.data.frame(Gentf)
}
Rf <- genotypefunc(Gen = "R")
dim(Rf)
## [1] 615 7
Kf <- genotypefunc(Gen = "K")
dim(Kf)
## [1] 501 7
Gf <- genotypefunc(Gen = "G")
dim(Gf)
## [1] 678 7
Yf <- genotypefunc(Gen = "Y")
dim(Yf)
## [1] 678 7
#combine together
GenRYGBK <- rbind(Rf, Kf, Bf, Gf, Yf)
dim(GenRYGBK)
## [1] 2970 7
Summarizing the data with genotypes
#merge with single mix
colnames(singmixf)[4] <- "count_all"
GenRYGBK_sm <- merge(singmixf, GenRYGBK, by = c("Date","Frame"))
#dimension check
dim(singmixf)
## [1] 624 4
dim(GenRYGBK)
## [1] 2970 7
dim(GenRYGBK_sm) #summarize - just B, K, R, Y, G
## [1] 2970 9
sumgen1 <- Rmisc::summarySE(GenRYGBK_sm, measurevar = c("prop"), groupvars = c("Date","Genotype","sing_mix", "Coral_stat"))
sum_onlygrybk <- Rmisc::summarySE(GenRYGBK_sm, measurevar = c("prop"), groupvars = c("Date","sing_mix", "Coral_stat"))
Model 2: Including genotypes in understanding disease prevalence
disease2 <- GenRYGBK_sm[which(GenRYGBK_sm$Coral_stat == "X"),]
disease2$Frame <- as.factor(disease2$Frame)
disease2$sing_mix <- as.factor(disease2$sing_mix)
disease2$Genotype <- as.factor(disease2$Genotype)
disease2$Date2 <- gsub("[.]", "/", disease2$Date)
disease2$Date3 <- as.Date(disease2$Date2)
#between genotype and single mixed
disease2_dens <- merge(Frame_area, disease2, by = "Frame")
dim(disease2)
## [1] 990 11
dim(disease2_dens)
## [1] 990 12
disease2_dens$density <- disease2_dens$countf/disease2_dens$Area
library(stats)
x2 <-poly(disease2_dens$Date3,2)
library("glmmTMB")
glmmtb_all <-
glmmTMB(cbind(count, countf-count) ~ sing_mix*x2*Genotype + density + (1|Frame), family = binomial, data= disease2_dens)
summary(glmmtb_all)
## Family: binomial ( logit )
## Formula:
## cbind(count, countf - count) ~ sing_mix * x2 * Genotype + density +
## (1 | Frame)
## Data: disease2_dens
##
## AIC BIC logLik deviance df.resid
## 2622.6 2779.4 -1279.3 2558.6 958
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Frame (Intercept) 1.611 1.269
## Number of obs: 990, groups: Frame, 26
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.01091 0.55771 -5.399 6.71e-08 ***
## sing_mixsingle 2.76577 1.14032 2.425 0.015290 *
## x21 70.04911 18.12525 3.865 0.000111 ***
## x22 -61.49875 18.21510 -3.376 0.000735 ***
## GenotypeG 0.09659 0.31877 0.303 0.761891
## GenotypeK 0.87877 0.37570 2.339 0.019335 *
## GenotypeR 0.62232 0.31792 1.957 0.050293 .
## GenotypeY 0.72232 0.32005 2.257 0.024013 *
## density -0.06856 0.12141 -0.565 0.572262
## sing_mixsingle:x21 -53.44289 18.56310 -2.879 0.003990 **
## sing_mixsingle:x22 27.44410 18.82341 1.458 0.144847
## sing_mixsingle:GenotypeG -2.37924 1.01933 -2.334 0.019589 *
## sing_mixsingle:GenotypeK -1.00748 1.04460 -0.964 0.334813
## sing_mixsingle:GenotypeR -1.16649 1.09148 -1.069 0.285193
## sing_mixsingle:GenotypeY -1.01770 1.11115 -0.916 0.359721
## x21:GenotypeG -42.77861 19.46394 -2.198 0.027961 *
## x22:GenotypeG 54.64573 19.68237 2.776 0.005497 **
## x21:GenotypeK 39.67457 24.72132 1.605 0.108522
## x22:GenotypeK -35.15077 25.40060 -1.384 0.166403
## x21:GenotypeR -79.63663 18.81900 -4.232 2.32e-05 ***
## x22:GenotypeR 33.73691 18.98273 1.777 0.075528 .
## x21:GenotypeY -54.68006 19.79114 -2.763 0.005730 **
## x22:GenotypeY 20.12562 20.57758 0.978 0.328056
## sing_mixsingle:x21:GenotypeG 42.88986 20.70042 2.072 0.038272 *
## sing_mixsingle:x22:GenotypeG -48.01809 21.37918 -2.246 0.024703 *
## sing_mixsingle:x21:GenotypeK 19.36679 25.44816 0.761 0.446640
## sing_mixsingle:x22:GenotypeK 6.29700 26.26124 0.240 0.810498
## sing_mixsingle:x21:GenotypeR 63.39830 19.74293 3.211 0.001322 **
## sing_mixsingle:x22:GenotypeR -50.79416 20.61550 -2.464 0.013744 *
## sing_mixsingle:x21:GenotypeY 63.16198 20.82326 3.033 0.002419 **
## sing_mixsingle:x22:GenotypeY -29.32972 21.98987 -1.334 0.182275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#anova(glm3a)
#car::Anova(glmmtb_all, type = "III")
Model 2 summary table
knitr::kable(car::Anova(glmmtb_all, type = "III"))
| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| (Intercept) | 29.1458553 | 1 | 0.0000001 |
| sing_mix | 5.8827825 | 1 | 0.0152896 |
| x2 | 15.0168983 | 2 | 0.0005484 |
| Genotype | 16.2412186 | 4 | 0.0027120 |
| density | 0.3189123 | 1 | 0.5722620 |
| sing_mix:x2 | 15.1217178 | 2 | 0.0005204 |
| sing_mix:Genotype | 5.8737111 | 4 | 0.2087802 |
| x2:Genotype | 109.9652765 | 8 | 0.0000000 |
| sing_mix:x2:Genotype | 39.1258261 | 8 | 0.0000047 |
library(DHARMa, quietly = TRUE)
glmtb_resid <- simulateResiduals(glmmtb_all)
testOutliers(glmtb_resid, type = "bootstrap")
##
## DHARMa bootstrapped outlier test
##
## data: glmtb_resid
## outliers at both margin(s) = 0, observations = 990, p-value = 0.72
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.008611111
## sample estimates:
## outlier frequency (expected: 0.00237373737373737 )
## 0
plot(glmtb_resid)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
#no significant outliers
Figure 2
disease2 <- disease2[which(disease2$Date!= "2019.02.01"),]
sum_gen_all_sing_mix2 <- Rmisc::summarySE(disease2, measurevar = c("prop"), groupvars = c("Date3","sing_mix"), na.rm = T)
#summary for Figure 2
sum_gybrk_sing_mix <- Rmisc::summarySE(na.omit(disease2), measurevar = c("prop"), groupvars = c("Date","Date3","sing_mix","Genotype"))
#summary for SFig 3
sumall_gybrk_sing_mix <- Rmisc::summarySE(na.omit(GenRYGBK_sm), measurevar = c("prop"), groupvars = c("Date","sing_mix","Genotype","Coral_stat"))
Supplement Figure S3: Proportion of health status by genotypes on both single and mixed frames
sumall_gybrk_sing_mix$Date2 <- gsub("[.]", "/", sumall_gybrk_sing_mix$Date)
sumall_gybrk_sing_mix$Date3 <- as.Date(sumall_gybrk_sing_mix$Date2)
sumall_gybrk_sing_mix$Coral_stat <- factor(sumall_gybrk_sing_mix$Coral_stat, levels = c("H","X","D"), labels = c("Healthy","Disease","Dead"))
ggplot(sumall_gybrk_sing_mix[which(sumall_gybrk_sing_mix$Date != "2019.02.01"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Proportion of Corals") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity")+ theme(text = element_text(size = 16), axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(vars(Genotype), vars(Coral_stat))
Individual genotype plots for disease prevalence
plot.Y <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="Y"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype) + lims(y = c(0,0.85))
plot.G <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="G"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))
plot.R <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="R"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))
plot.B <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="B"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))
plot.K <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="K"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))
Figure 2
library(patchwork)
(sing_mix_all + plot_layout(widths = 25))/(plot.G + plot.Y + plot.R + plot.B + plot.K
+ plot_layout(nrow = 1, widths = 25))/plot_layout(guides = "collect") & theme(legend.position = "bottom") & plot_annotation(tag_levels = 'a')
Supplemental Figure 1: Temperature
tempcsv <- read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Microbial Sampling_wholeframes/Temperature/tempcsv.csv", header=FALSE, comment.char="#")
temp <- tempcsv[,2:3]
colnames(temp) <- c("Date_time","temp")
temp <- separate(temp, "Date_time",c("Date","time")," ")
temp$Date <- as.Date(temp$Date, format = "%m/%d/%y")
temp2 <- subset(temp, Date < "2019-08-19")
temp.sum <- summarySE(temp2, measurevar = "temp", groupvars = c("Date"))
martha_r <- read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Microbial Sampling_wholeframes/Temperature/martha_r.csv")
temp3 <-martha_r
colnames(temp3) <- c("Date_time","temp")
temp4 <- separate(temp3, "Date_time",c("Date","time")," ")
str(temp4)
## 'data.frame': 11179 obs. of 3 variables:
## $ Date: chr "2/11/19" "2/11/19" "2/11/19" "2/11/19" ...
## $ time: chr "15:30" "16:00" "16:30" "17:00" ...
## $ temp: num 80.7 80.7 80.6 80.6 80.6 ...
temp4$Date <- as.Date(temp4$Date, format = "%m/%d/%y")
temp5 <- subset(temp4, Date > "2019-05-19" & Date < "2019-08-10")
temp.sum2 <- summarySE(temp5, measurevar = "temp", groupvars = c("Date"))
dim(temp.sum2)
## [1] 82 6
temp.sum2$Loc <- "Martha"
temp.sum$Loc <- "Deep"
temp.both <- data.table(rbind(temp.sum, temp.sum2))
temp.both$tempC <- (temp.both$temp - 32)* (5/9)
templab <- "Temperature (°C)"
templot <- ggplot(temp.both, aes(x = Date, y = tempC)) + geom_point(size = 3, fill = "black", color = "black", pch = 21) + scale_x_date(date_labels = "%m-%d",date_breaks = "10 day" ) + ylab(templab) + theme(text = element_text(size = 16)) + theme_bw() + geom_line()
library(cowplot)
sum_all2_a <- select(sum_all2, Coral_stat, N, prop, se, Date3)
colnames(sum_all2_a)[5] <- "Date"
temp.dis <- merge(sum_all2_a, temp.both, by = "Date")
tempdeadp <- ggplot(temp.dis[which(temp.dis$Coral_stat == "D"),], aes(x = tempC, y = prop)) + geom_point(size = 3) + theme_bw() + xlab(templab) + ylab("Proportion of Dead Corals") + geom_smooth(se = F, col = "black")
tempdisp <-ggplot(temp.dis[which(temp.dis$Coral_stat == "X"),], aes(x = tempC, y = prop)) + geom_point(size = 3) + theme_bw() + xlab(templab) + ylab("Proportion of Diseased Corals") + geom_smooth(se = F, col = "gray")
healthtemp <- ggplot(temp.dis, aes(x = tempC, y = prop, group= Coral_stat)) + geom_point(size = 3, aes(color = Coral_stat)) + theme_bw() + xlab(templab) + ylab("Proportion of Corals") + geom_smooth(se = F, aes(color = Coral_stat), method = "lm")+ scale_color_manual(values = c("black","orange", "gray"), name = "Health", labels = c("Dead","Healthy", "Disease"))
plot_grid(templot, healthtemp, labels = "auto")
## `geom_smooth()` using formula 'y ~ x'
temp.prop.lm <- lm(prop~tempC*Coral_stat, data= temp.dis)
summary(temp.prop.lm)
##
## Call:
## lm(formula = prop ~ tempC * Coral_stat, data = temp.dis)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.055072 -0.013532 0.002014 0.013404 0.053653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.936543 0.266787 -7.259 1.57e-09 ***
## tempC 0.068640 0.008957 7.664 3.45e-10 ***
## Coral_statH 5.796832 0.377293 15.364 < 2e-16 ***
## Coral_statX 1.012797 0.377293 2.684 0.00963 **
## tempC:Coral_statH -0.179529 0.012667 -14.173 < 2e-16 ***
## tempC:Coral_statX -0.026392 0.012667 -2.084 0.04195 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02192 on 54 degrees of freedom
## Multiple R-squared: 0.988, Adjusted R-squared: 0.9869
## F-statistic: 891 on 5 and 54 DF, p-value: < 2.2e-16
plot(temp.prop.lm)
car::Anova(temp.prop.lm, type = "III")
## Anova Table (Type III tests)
##
## Response: prop
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.025308 1 52.69 1.565e-09 ***
## tempC 0.028209 1 58.73 3.452e-10 ***
## Coral_stat 0.129381 2 134.68 < 2.2e-16 ***
## tempC:Coral_stat 0.112518 2 117.13 < 2.2e-16 ***
## Residuals 0.025937 54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(temp.prop.lm)
## (Intercept) tempC Coral_statH Coral_statX
## -1.93654313 0.06864041 5.79683246 1.01279693
## tempC:Coral_statH tempC:Coral_statX
## -0.17952941 -0.02639181
temp.prop.lmX <- lm(prop~tempC, data= temp.dis[which(temp.dis$Coral_stat == "X"),])
temp.prop.lmD <- lm(prop~tempC, data= temp.dis[which(temp.dis$Coral_stat == "D"),])
temp.prop.lmH <- lm(prop~tempC, data= temp.dis[which(temp.dis$Coral_stat == "H"),])
summary(temp.prop.lmX)
##
## Call:
## lm(formula = prop ~ tempC, data = temp.dis[which(temp.dis$Coral_stat ==
## "X"), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.055072 -0.013532 0.004643 0.018411 0.037225
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.92375 0.30080 -3.071 0.006583 **
## tempC 0.04225 0.01010 4.184 0.000558 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02471 on 18 degrees of freedom
## Multiple R-squared: 0.493, Adjusted R-squared: 0.4648
## F-statistic: 17.5 on 1 and 18 DF, p-value: 0.0005581
summary(temp.prop.lmD)
##
## Call:
## lm(formula = prop ~ tempC, data = temp.dis[which(temp.dis$Coral_stat ==
## "D"), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0177993 -0.0058061 0.0007323 0.0074983 0.0166745
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.936543 0.130610 -14.83 1.57e-11 ***
## tempC 0.068640 0.004385 15.65 6.30e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01073 on 18 degrees of freedom
## Multiple R-squared: 0.9316, Adjusted R-squared: 0.9278
## F-statistic: 245 on 1 and 18 DF, p-value: 6.301e-12
summary(temp.prop.lmH)
##
## Call:
## lm(formula = prop ~ tempC, data = temp.dis[which(temp.dis$Coral_stat ==
## "H"), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04191 -0.01733 0.00134 0.01598 0.05365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.86029 0.32556 11.86 6.12e-10 ***
## tempC -0.11089 0.01093 -10.14 7.15e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02674 on 18 degrees of freedom
## Multiple R-squared: 0.8512, Adjusted R-squared: 0.8429
## F-statistic: 102.9 on 1 and 18 DF, p-value: 7.148e-09